home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / multiroots / test.c < prev   
Encoding:
C/C++ Source or Header  |  2001-08-22  |  7.1 KB  |  261 lines
  1. /* multiroots/test.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <stdio.h>
  23. #include <gsl/gsl_vector.h>
  24. #include <gsl/gsl_test.h>
  25. #include <gsl/gsl_multiroots.h>
  26.  
  27. #include <gsl/gsl_ieee_utils.h>
  28.  
  29. #include "test_funcs.h"
  30. int test_fdf (const char * desc, gsl_multiroot_function_fdf * function, initpt_function initpt, double factor, const gsl_multiroot_fdfsolver_type * T);
  31. int test_f (const char * desc, gsl_multiroot_function_fdf * fdf, initpt_function initpt, double factor, const gsl_multiroot_fsolver_type * T);
  32.  
  33.  
  34. int 
  35. main (void)
  36. {
  37.   const gsl_multiroot_fsolver_type * fsolvers[5] ;
  38.   const gsl_multiroot_fsolver_type ** T1 ;
  39.  
  40.   const gsl_multiroot_fdfsolver_type * fdfsolvers[5] ;
  41.   const gsl_multiroot_fdfsolver_type ** T2 ;
  42.  
  43.   double f;
  44.  
  45.   fsolvers[0] = gsl_multiroot_fsolver_dnewton;
  46.   fsolvers[1] = gsl_multiroot_fsolver_broyden;
  47.   fsolvers[2] = gsl_multiroot_fsolver_hybrid;
  48.   fsolvers[3] = gsl_multiroot_fsolver_hybrids;
  49.   fsolvers[4] = 0;
  50.  
  51.   fdfsolvers[0] = gsl_multiroot_fdfsolver_newton;
  52.   fdfsolvers[1] = gsl_multiroot_fdfsolver_gnewton;
  53.   fdfsolvers[2] = gsl_multiroot_fdfsolver_hybridj;
  54.   fdfsolvers[3] = gsl_multiroot_fdfsolver_hybridsj;
  55.   fdfsolvers[4] = 0;
  56.  
  57.   gsl_ieee_env_setup();
  58.  
  59.  
  60.   f = 1.0 ;
  61.   
  62.   T1 = fsolvers ;
  63.   
  64.   while (*T1 != 0) 
  65.     {
  66.       test_f ("Rosenbrock", &rosenbrock, rosenbrock_initpt, f, *T1);
  67.       test_f ("Roth", &roth, roth_initpt, f, *T1);
  68.       test_f ("Powell badly scaled", &powellscal, powellscal_initpt, f, *T1);
  69.       test_f ("Brown badly scaled", &brownscal, brownscal_initpt, f, *T1);
  70.       test_f ("Powell singular", &powellsing, powellsing_initpt, f, *T1);
  71.       test_f ("Wood", &wood, wood_initpt, f, *T1);
  72.       test_f ("Helical", &helical, helical_initpt, f, *T1);
  73.       test_f ("Discrete BVP", &dbv, dbv_initpt, f, *T1);
  74.       test_f ("Trig", &trig, trig_initpt, f, *T1);
  75.       T1++;
  76.     }
  77.   
  78.   T2 = fdfsolvers ;
  79.   
  80.   while (*T2 != 0) 
  81.     {
  82.       test_fdf ("Rosenbrock", &rosenbrock, rosenbrock_initpt, f, *T2);
  83.       test_fdf ("Roth", &roth, roth_initpt, f, *T2);
  84.       test_fdf ("Powell badly scaled", &powellscal, powellscal_initpt, f, *T2);
  85.       test_fdf ("Brown badly scaled", &brownscal, brownscal_initpt, f, *T2);
  86.       test_fdf ("Powell singular", &powellsing, powellsing_initpt, f, *T2);
  87.       test_fdf ("Wood", &wood, wood_initpt, f, *T2);
  88.       test_fdf ("Helical", &helical, helical_initpt, f, *T2);
  89.       test_fdf ("Discrete BVP", &dbv, dbv_initpt, f, *T2);
  90.       test_fdf ("Trig", &trig, trig_initpt, f, *T2);
  91.       T2++;
  92.     }
  93.  
  94.   exit (gsl_test_summary ());
  95. }
  96.  
  97. void scale (gsl_vector * x, double factor);
  98.  
  99. void
  100. scale (gsl_vector * x, double factor)
  101. {
  102.   size_t i, n = x->size;
  103.  
  104.   if (gsl_vector_isnull(x))
  105.     {
  106.       for (i = 0; i < n; i++)
  107.         {
  108.           gsl_vector_set (x, i, factor);
  109.         }
  110.     }
  111.   else
  112.     {
  113.       for (i = 0; i < n; i++)
  114.         {
  115.           double xi = gsl_vector_get(x, i);
  116.           gsl_vector_set(x, i, factor * xi);
  117.         }
  118.     } 
  119. }
  120.  
  121. int
  122. test_fdf (const char * desc, gsl_multiroot_function_fdf * function, 
  123.           initpt_function initpt, double factor,
  124.           const gsl_multiroot_fdfsolver_type * T)
  125. {
  126.   int status;
  127.   double residual = 0;
  128.   size_t i, n = function->n, iter = 0;
  129.   
  130.   gsl_vector *x = gsl_vector_alloc (n);
  131.   gsl_matrix *J = gsl_matrix_alloc (n, n);
  132.  
  133.   gsl_multiroot_fdfsolver *s;
  134.  
  135.   (*initpt) (x);
  136.  
  137.   if (factor != 1.0) scale(x, factor);
  138.  
  139.   s = gsl_multiroot_fdfsolver_alloc (T, n);
  140.   gsl_multiroot_fdfsolver_set (s, function, x);
  141.  
  142.   do
  143.     {
  144.       iter++;
  145.       status = gsl_multiroot_fdfsolver_iterate (s);
  146.       
  147.       if (status)
  148.         break ;
  149.  
  150.       status = gsl_multiroot_test_residual (s->f, 0.0000001);
  151.     }
  152.   while (status == GSL_CONTINUE && iter < 1000);
  153.  
  154. #ifdef DEBUG
  155.   printf("x "); gsl_vector_fprintf (stdout, s->x, "%g"); printf("\n");
  156.   printf("f "); gsl_vector_fprintf (stdout, s->f, "%g"); printf("\n");
  157. #endif
  158.  
  159.  
  160. #ifdef TEST_JACOBIAN
  161.  {
  162.     double r,sum; size_t j;
  163.  
  164.     gsl_multiroot_function f1 ;
  165.     f1.f = function->f ;
  166.     f1.n = function->n ;
  167.     f1.params = function->params ;
  168.     
  169.     gsl_multiroot_fdjacobian (&f1, s->x, s->f, GSL_SQRT_DBL_EPSILON, J);
  170.   
  171.     /* compare J and s->J */
  172.     
  173.     r=0;sum=0;
  174.     for (i = 0; i < n; i++)
  175.       for (j = 0; j< n ; j++)
  176.         {
  177.           double u = gsl_matrix_get(J,i,j);
  178.           double su = gsl_matrix_get(s->J, i, j);
  179.           r = fabs(u - su)/(1e-6 + 1e-6 * fabs(u)); sum+=r;
  180.           if (fabs(u - su) > 1e-6 + 1e-6 * fabs(u))
  181.             printf("broken jacobian %g\n", r);
  182.         }
  183.     printf("avg r = %g\n", sum/(n*n));
  184.   }
  185. #endif
  186.  
  187.   for (i = 0; i < n ; i++)
  188.     {
  189.       residual += fabs(gsl_vector_get(s->f, i));
  190.     }
  191.  
  192.   gsl_multiroot_fdfsolver_free (s);
  193.   gsl_matrix_free(J);
  194.   gsl_vector_free(x);
  195.  
  196.   gsl_test(status, "%s on %s (%g), %u iterations, residual = %.2g", T->name, desc, factor, iter, residual);
  197.  
  198.   return status;
  199. }
  200.  
  201.  
  202. int
  203. test_f (const char * desc, gsl_multiroot_function_fdf * fdf, 
  204.         initpt_function initpt, double factor,
  205.         const gsl_multiroot_fsolver_type * T)
  206. {
  207.   int status;
  208.   size_t i, n = fdf->n, iter = 0;
  209.   double residual = 0;
  210.  
  211.   gsl_vector *x;
  212.  
  213.   gsl_multiroot_fsolver *s;
  214.   gsl_multiroot_function function;
  215.  
  216.   function.f = fdf->f;
  217.   function.params = fdf->params;
  218.   function.n = n ;
  219.  
  220.   x = gsl_vector_alloc (n);
  221.  
  222.   (*initpt) (x);
  223.  
  224.   if (factor != 1.0) scale(x, factor);
  225.  
  226.   s = gsl_multiroot_fsolver_alloc (T, n);
  227.   gsl_multiroot_fsolver_set (s, &function, x);
  228.  
  229. /*   printf("x "); gsl_vector_fprintf (stdout, s->x, "%g"); printf("\n"); */
  230. /*   printf("f "); gsl_vector_fprintf (stdout, s->f, "%g"); printf("\n"); */
  231.  
  232.   do
  233.     {
  234.       iter++;
  235.       status = gsl_multiroot_fsolver_iterate (s);
  236.       
  237.       if (status)
  238.         break ;
  239.  
  240.       status = gsl_multiroot_test_residual (s->f, 0.0000001);
  241.     }
  242.   while (status == GSL_CONTINUE && iter < 1000);
  243.  
  244. #ifdef DEBUG
  245.   printf("x "); gsl_vector_fprintf (stdout, s->x, "%g"); printf("\n");
  246.   printf("f "); gsl_vector_fprintf (stdout, s->f, "%g"); printf("\n");
  247. #endif
  248.  
  249.   for (i = 0; i < n ; i++)
  250.     {
  251.       residual += fabs(gsl_vector_get(s->f, i));
  252.     }
  253.  
  254.   gsl_multiroot_fsolver_free (s);
  255.   gsl_vector_free(x);
  256.  
  257.   gsl_test(status, "%s on %s (%g), %u iterations, residual = %.2g", T->name, desc, factor, iter, residual);
  258.  
  259.   return status;
  260. }
  261.